
%% runParamSingle.m 
% 
% This code loads a mode from an XLS-07 file, solves the model and displays
% a number of theoretical moments and plots the IRFS for the observables. 
% Use to debug against Dynare code 
% 
% See *runPeramVariation*.m for a code that does this the same for a
% perturbation of the parameters 
%
% Alejandro Justiniano (c) Oct 30 2013
%=========================================================================
clear all; 
close all; 
clc; 

cucd=cd; 

% =========================================================================
%% INPUTS
% Function Handle 
funcHandle=@modBMJ_newBaseDepUNomShare; 
addsol.GDPCWCons=1; 
addsol.measureEQ.wage.numSeries     = 2;
addsol.measureEQ.wage.nameSeries    = {'LEPRIVA_Core' 'LS_Core'};
addsol.measureEQ.price.numSeries    = 3;		
addsol.measureEQ.price.nameSeries   = {'JCXFEBM_LD100'	'JCMXFEM_LD100' 'PCUSLFEN_LD100SA'};

% Name and location of the XLS-07 filename that contains the mode 
xls.filename='priors_firstmom 11 25 14.xls'; 
xls.sheet='ModePrecision'; 
xls.path=fullfile(cucd,'BMJ_new','BaseDepUFirstMom','FirstMoments_S1','priors_firstmom 11 25 14'); 

% Number of horizons for the IRFS 
dim.horIRF=12*40; 
% Number of horizons for the forecast error decompositions 
dim.horDec=20; 

%==========================================================================
%% Load XLS File with Mode 
cd(xls.path);
if isempty(xls.sheet)==false
    [parVec,parNames]=xlsread(xls.filename,xls.sheet);
else
    [parVec,parNames]=xlsread(xls.filename,xls.sheet);
end
cd(cucd); 

%==========================================================================    
%% Solve Model 
[GG,RR,CONS,eu,SDX,ZZ,~,~,flags,~,stateNames,shockNames]=feval(funcHandle,parVec,[],addsol); 
if isstruct(shockNames)==true 
    shockNames=shockNames.short; 
end 
%% Extract the names of the observables as defined in ZZ 
[obsNames,posInStates]=extractObsNames(ZZ(:,:,end),stateNames);  

% Observables
dim.y=size(ZZ,1); 
% Full state of transition equation 
dim.s=size(GG,1); 
% Innovations 
dim.eta=size(SDX,1); 

%% Compute all Theoretical, i.e. asymptotic moments 
% uses old code...pre-structure architechture
[stdObs,stdStates,irfObs,irfStates,ACZeroObs,ACZOneObs,vDecObs,vDecStates,vDecHorObs,vDecHorStates]=...
    irfandthmom(GG,RR,SDX,ZZ,dim.horIRF,dim.horDec,zeros(size(ZZ,1),1),0); 

%% 4, Simple code to plot impulse responses 
%  IRF stored in irfObs [dim.y dim.horIRF dim.eta] 
%  i.e, irfObs(ii,:,jj) is IRF of observable series ii to shock jj 
%  Magnitude is 1 std. dev. impulse 
%{
fig.vec=zeros(dim.eta,1); 
fig.cols=3; 
fig.rows=ceil( dim.y / fig.cols ); 
fig.xAxis=(0:dim.horIRF); 
fig.xAxis=fig.xAxis(:); 

for shind=1:dim.eta 
    
    fig.vec(shind)=figure; 
    
    for serind=1:dim.y 
        
        subplot( fig.rows, fig.cols, serind ); 
        
        plot(fig.xAxis,squeeze( irfObs(serind,:,shind) ),'b','LineWidth',1.5); 
        
        title(char(obsNames{serind})); 
        
        set(gca,'box','off'); 
        
        axis tight; 
        
    end

    set(fig.vec(shind),'PaperPosition',[0.75 0.5 7.5 10]);

    suptitle(['Shock: ',char(shockNames{shind})],14);
end         
%}
        
%% 4, Simple code to plot impulse responses (States)
%  IRF stored in irfObs [dim.y dim.horIRF dim.eta] 
%  i.e, irfObs(ii,:,jj) is IRF of observable series ii to shock jj 
%  Magnitude is 1 std. dev. impulse 
PlotStates = {'Consumption' 'GDP' 'k physical (k)' 'r^k'};
PlotStatesLoc = cellposition(PlotStates,stateNames)';

PlotShocks = {'Government Drift'};
PlotShocksLoc = cellposition(PlotShocks,shockNames)';

fig.vec=zeros(numel(PlotShocks),1); 
fig.cols=1; 
fig.rows=4; 
fig.xAxis=(0:dim.horIRF); 
fig.xAxis=fig.xAxis(:); 


figureind = 1;
for shind=PlotShocksLoc
    
    fig.vec(figureind)=figure; 
    
    subplotind = 1;
    for serind=PlotStatesLoc
        
        subplot(fig.rows,fig.cols,subplotind); 
        
        plot(fig.xAxis,-1*squeeze(irfStates(serind,:,shind)),'b','LineWidth',1.5); 
        
        title(char(stateNames{serind})); 
        
        set(gca,'box','off'); 
        
        axis tight; 
        subplotind = subplotind+1;
        
    end

    set(fig.vec(figureind),'PaperPosition',[0.75 0.5 7.5 10]);

    suptitle(['Shock: ',char(shockNames{shind})],14);

    figureind = figureind+1;
end         
      